#' ---
#' title: "Perceived Inequality and Populism"
#' description: "Replication Observational Evidence"
#' author: "Lukas F. Stoetzer"
#' date: "June 2025"
#' ---

# load packages
library(tidyverse)
library(sandwich)
library(lmtest)
library(MatchIt)
library(texreg)

# Load Cleaned Data
issp <- readRDS("00_Data_ISSP.RDS")


# Estimation function ==============

  # Prepare Data
  issp_any <- issp %>% 
    filter(POPPARTY == "ANY") %>%
    mutate(AGE_CAT = cut(AGE,breaks = 5),
           CNTRY_ISSP = paste(CNTRY,ISSP,sep="_")) %>% 
    select("SEX", "AGE_CAT","AGE", "INC_CAT", "ISSP",
           "DEGREE_UNI","DEGREE", "UNEMPLOYED","WEIGHT",
           "CNTRY_ISSP", "CNTRY", "ISSP", INQ, INQ_PYR, VOTE) %>% 
    na.exclude()

  # Equations (with controls)  
  eqn1 <- paste("VOTE ~ INQ + CNTRY_ISSP +  
                  SEX + AGE_CAT + INC_CAT + DEGREE_UNI + UNEMPLOYED") 
  eqn2 <- paste("VOTE ~ as.factor(INQ_PYR) + CNTRY_ISSP + 
                  SEX + AGE_CAT + INC_CAT + DEGREE_UNI + UNEMPLOYED") 
  
  # Linear Models
  res.lm1 <- lm(eqn1,issp_any)
  res.lm2 <- lm(eqn2,issp_any)
  
  # Matching
  match_out <- matchit(INQ ~ 
                         CNTRY_ISSP + SEX + 
                         AGE_CAT + INC_CAT + DEGREE_UNI + UNEMPLOYED, 
                       data=issp_any, method = "exact")
  
  # Matched data
  dat_match <- match.data(match_out)
  
  # Balance Table
  # Extract balance summary
  balance_summary <- summary(match_out, un = TRUE)
  
  # Get balance data before and after matching
  balance_data <- data.frame(
    Means_Before_Control = balance_summary$sum.all[,"Means Treated"]*100,
    Means_Before_Treated = balance_summary$sum.all[,"Means Control"]*100,
    Means_Before_Diff = balance_summary$sum.all[,"Std. Mean Diff."]*100, 
    
    Means_After_Control = balance_summary$sum.matched[,"Means Treated"]*100,
    Means_After_Treated = balance_summary$sum.matched[,"Means Control"]*100,
    Means_After_Diff = balance_summary$sum.matched[,"Std. Mean Diff."]*100
  )
  
  # Models
  res.exact1 <- lm(eqn1,data = dat_match, weights = weights)
  res.exact2 <- lm(eqn2,data = dat_match, weights = weights)
  
  # Put in list
  mdls_pooled <- list(res.lm1,res.lm2,res.exact1,res.exact2)

  # Regression Table
  htmlreg(mdls_pooled,digits = 3,ci.force = T, 
          file = "Table1.html",
            custom.coef.map = list(
              "INQ" = "Preceived Inequality Type A",
              "as.factor(INQ_PYR)2" = "... Type B",
              "as.factor(INQ_PYR)3" = "... Type C",
              "as.factor(INQ_PYR)4" = "... Type D",
              "as.factor(INQ_PYR)5" = "... Type E"
            ))

  
# Model for parties separately ===========

  # Prepare Data
  issp_party <- issp %>% 
    filter(POPPARTY != "ANY") %>%
    mutate(CASE = paste0(CNTRY,"_",POPPARTY,"_",ISSP,sep="")) 
  
  # Unique parties
  unique_parties <- issp_party %>%
    pull(CASE) %>% unique()

  # Loop over unique party-country cases
  res_party <- list()
  for(prty in  unique_parties){
    
    # filter data
    dat <-  issp_party %>% 
      filter(CASE == prty) %>% 
      mutate(AGE_CAT = cut(AGE,breaks = 5),
             CNTRY_ISSP = paste(CNTRY,ISSP,sep="_")) %>% 
      select("SEX", "AGE_CAT","AGE", "INC_CAT", "ISSP",
             "DEGREE_UNI","DEGREE", "UNEMPLOYED","WEIGHT",
             "CNTRY_ISSP", "CNTRY", "ISSP", INQ, INQ_PYR, VOTE) %>% 
      na.exclude()
      
    # Equations (with controls)  
    eqn1 <- paste("VOTE ~ INQ +   
                  SEX + AGE_CAT + INC_CAT + DEGREE_UNI + UNEMPLOYED") 
    eqn2 <- paste("VOTE ~ as.factor(INQ_PYR) +  
                  SEX + AGE_CAT + INC_CAT + DEGREE_UNI + UNEMPLOYED") 
    
    # Linear Models
    res.lm1 <- lm(eqn1,dat)
    res.lm2 <- lm(eqn2,dat)
    
    # Matching
    match_out <- matchit(INQ ~ SEX + AGE_CAT + INC_CAT + DEGREE_UNI + UNEMPLOYED, 
                         data=dat, method = "exact")
    
    # Matched data
    dat_match <- match.data(match_out)
    
    # Models
    res.exact1 <- lm(eqn1,data = dat_match, weights = weights)
    res.exact2 <- lm(eqn2,data = dat_match, weights = weights)
    
    # save in list 
    res_party[[prty]] <-  list("lm.m1" = res.lm1, "lm.m2"=res.lm2,
                               "exact.m1" = res.exact1, "exact.m2" =res.exact2)
      
  }
  
  
  # Get estimates
  get_est <-  function(i, what="lm.m1") {
    
    
    # Get names
    var_nam <- names(coef(res_party[[i]][[what]]))
    var_sel <- grep("INQ",var_nam)
    
    
    # Create Data Frame
    df <- data.frame( 
      "mdl" = what,
      "var_nam" = var_nam[var_sel],
      "case" = names(res_party)[[i]],
      "est"=coef(res_party[[i]][[what]])[var_sel],
      "se" =sqrt(diag(vcov(res_party[[i]][[what]]))[var_sel])) %>%
      mutate(
        est_low = est - 1.96*se,
        est_hgh = est + 1.96*se,
        sig = (est/se > 1.96)
      )
    
    rownames(df) <- NULL
    
    return(df)
  }
  
  # For models
  res_est_party <- bind_rows(
  lapply(seq_along(res_party),get_est,what="lm.m1"),
  lapply(seq_along(res_party),get_est,what="lm.m2"),
  lapply(seq_along(res_party),get_est,what="exact.m1"),
  lapply(seq_along(res_party),get_est,what="exact.m2")) %>%
    separate(case, c("CNTRY","PARTY","ISSP")) 
  
  # Add left-right info
  res_est_party <- res_est_party %>%
    mutate(PARTY = ifelse(PARTY == "SNSSL", "SNS", PARTY)) %>%
    left_join(.,select(read_csv("00_Data_Parties.csv"),-ISSP)) %>%
    mutate(CntryParty = paste(PARTY, " (", CNTRY, " " ,ISSP, ")", sep=""),
         LEFT_RIGHT_DEU = as.factor(case_when(
           LEFT_RIGHT == "FarRight" ~ "Far-Right Populist \n Party",
           LEFT_RIGHT == "FarLeft" ~ "Far-Left Populist \n Party",
           LEFT_RIGHT == "None" ~ "Populist \n Party"
         )),
         LEFT_RIGHT_DEU = relevel(LEFT_RIGHT_DEU, "Far-Left Populist \n Party"))

  # Figure Model lm.m1
  ggplot(filter(res_est_party, mdl=="lm.m1")) +
    geom_pointrange(mapping = aes(x = reorder(CntryParty,est), y = est, ymin = est_low, ymax = est_hgh, col=!sig)) +
    scale_color_grey() +
    coord_flip() + 
    ylim(-0.3,0.4) +
    geom_hline(yintercept = 0, color = "red",alpha=0.5) +
    facet_wrap(~ LEFT_RIGHT_DEU, scales="free") +
    theme_minimal() + theme(legend.position = "none") +
    labs(y = "", x = "")

  ggsave("Figure2.pdf", width=9,height = 5)

  